Workflow: Capital Asset Pricing Model for Ecological Portfolios
Mark D. Scheuerell
Fish Ecology Division, Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA USA, mark.scheuerell@noaa.gov
0.0.1 Version
This is version 0.17.08.01.
1 Initialization
library(readr)
library(reshape2)
library(MARSS)2 Data
For these analyses we will use time series of counts of adult spring/summer Chinook salmon from the Snake River basin in Oregon and Idaho. This group of fish forms an Evolutionarily Significant Unit, which was listed as threatened under the US Endangered Species Act in 1992.
The data were compiled largely through the efforts of Tom Cooney (Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Portland, OR USA, tom.cooney@noaa.gov).
2.1 Assets
The assets data file has the following columns of information:
pop= abbreviated name of the populationcode= 5-letter abbreviation for MPG (2 letters) & population (3 letters)mpg= name of 1 of 5 Major Population Groupsha= estimated area (hectares) of available spawning habitatcategory= hatchery supplementation category (“impact” or “control”)brood_yr= year when spawning occurred
nS= total number of spawning fishphos= proportion of hatchery-origin spawners
rawdata <- read.csv("https://raw.githubusercontent.com/mdscheuerell/CAPM/master/data/srss_chin.csv")
head(rawdata)
## pop code mpg ha category brood_yr nS phos
## 1 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1955 122 0
## 2 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1956 2606 0
## 3 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1957 1754 0
## 4 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1958 486 0
## 5 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1959 712 0
## 6 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1960 3161 0The first thing to do is remove hatchery-origin spawners from the counts, and then convert the counts to log-density based on the estimated spawning area for each population. In a couple of years, incomplete censuses with small sample sizes led to estimated fractions of wild fish equal to zero. I will treat those as NA.
dat <- rawdata[,c("mpg", "code", "pop", "category", "brood_yr")]
## calculate log-density
dat$Sw <- round(log(rawdata$nS*(1-rawdata$phos)/rawdata$ha), 2)
## replace possible -Inf with NA
dat$Sw[is.infinite(dat$Sw)] <- NAFor easier bookkeeping, I’ll also switch the Imnaha’s MPG over to the Grande Ronde.
## change Imnaha name for easier sorting later
dat[,"code"] <- sub("IRMAI", "GR_IR", dat[,"code"], fixed=TRUE)
## MPG abbreviations
mpg_ID <- c("GR", "MF", "SF", "SR")Only a few of the populations have observations going all the way back to the 1950s, so I’ll trim the data accordingly.
## use generation time of 4 years
gt <- 4
## first brood year to consider
yr_first <- 1960
## last brood year to consider
yr_last <- 2011
## years of interest
t_index <- seq(yr_first, yr_last-gt)
## length of time period
TT <- length(t_index)
## subset data
dat <- subset(dat, brood_yr>=yr_first & brood_yr<=yr_last)I will use the MARSS package to fit the CAPM model, which requires the data to be formatted as an [N x T] matrix.
## reshape assets
assets <- dcast(dat, brood_yr ~ code, value.var="Sw")[,-1]
## compute differences
assets <- t(apply(assets, 2, diff, lag=gt))
matplot(t(assets), type="b")
## number of popns
NN <- dim(assets)[1]
## names of popns
names_pop <- rownames(assets)Finally, here are some nicer, longer names for the various populations.
# set up naming & ordering scheme
longnames <- c(
"Wenaha",
"Grande Ronde",
"Catherine",
"Minam",
"Lostine",
"Imnaha",
"SF Salmon",
"Secesh",
"EFSF Salmon",
"Big",
"Sulphur",
"Bear Valley",
"Marsh",
"Loon",
"Camas",
"Yankee Fork",
"Valley",
"UM Salmon",
"LM Salmon",
"EF Salmon",
"Lemhi")
name_ord <- c(6,3,5,4,2,1,12,10,15,14,13,11,9,7,8,20,21,19,18,17,16)
names_mat <- data.frame(ord=name_ord, ID=names_pop)
names_mat <- names_mat[order(names_mat$ord),]
names_mat$long <- longnames2.2 Market index
The first “market”" index is based on the monthly North Pacific Gyre Oscillation (NPGO) data provided by Emanuele Di Lorenzo, which are available here. We begin by downloading the raw NPGO data and viewing the metadata.
## URL for NPGO data
url_NPGO <- "http://www.o3d.org/npgo/npgo.php"
## raw NPGO data
all_NPGO <- read_lines(url_NPGO)
## line with data headers
hdr_NPGO <- which(lapply(all_NPGO,grep,pattern="YEAR")==1, arr.ind=TRUE)
## print PDO metadata
print(all_NPGO[seq(hdr_NPGO)],quote=FALSE)## [1]
## [2] <html>
## [3] <body>
## [4]
## [5] <pre># Last update 29-Mar-2017 by E. Di Lorenzo
## [6] # NPGO index monthly averages
## [7] # from Jan-1950 to Feb-2017
## [8] #
## [9] # WARNING: Values after Dec-2004 are updated
## [10] # using Satellite SSHa from AVISO Delayed Time product.
## [11] # http://opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-dt-global-allsat-msla-h
## [12] #
## [13] # PRELIMINARY: Values after Jan-2016 are preliminary and updated
## [14] # using Satellite SSHa from AVISO Near Real Time product.
## [15] # http://opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-nrt-over30d-global-allsat-msla-h
## [16] #
## [17] # The update is performed by taking the NPGO spatial pattern of Di Lorenzo et al. 2008
## [18] # computed over the period 1950-2004, and projecting the AVISO Satellite SSHa.
## [19] # During the pre-processing of the AVISO data, we remove the seasonal cycle based on
## [20] # the 1993-2004 seasonal means.
## [21] #
## [22] # AVISO PRODUCT UPDATE Summer 2014: AVISO has released a re-processed dataset for the sea level.
## [23] # Starting from the November 2014, the NPGO index is computed with this updated dataset. NPGO
## [24] # values from 2004 onward have been recomputed with very minor differences from previous releases.
## [25] #
## [26] # Ref:
## [27] # Di Lorenzo et al., 2008: North Pacific Gyre Oscillation
## [28] # links ocean climate and ecosystem change, GRL.
## [29] #
## [30] # YEAR MONTH NPGO indexNext, we will extract the actual NPGO indices for the years of interest and inspect the file contents. We also want the average annual NPGO index from January 1 through December 31 during the first year that the juvenile salmon are in the ocean (i.e., during their second year of life). Therefore, we need NPGO values from yr_first + 2 = 1962 through yr_last + 2 = 2013.
## NPGO data for years of interest
raw_NPGO <- read_table(url_NPGO, col_names = FALSE, skip=hdr_NPGO + (yr_first+2-1950)*12, n_max = TT*12)
colnames(raw_NPGO) <- c("year","month","NPGO")
## calculate annual means
dat_NPGO <- aggregate(raw_NPGO$NPGO, by = list(year = raw_NPGO$year), FUN = mean)
colnames(dat_NPGO) <- c("brood_yr", "NPGO")
## match calendar yr to brood year
dat_NPGO$brood_yr <- dat_NPGO$brood_yr - 2
dat_NPGO
## market index
market <- round(dat_NPGO$NPGO, 3)
plot.ts(market)
cor(market, t(assets), use="pairwise.complete.obs")## brood_yr NPGO
## 1 1960 -0.093424126
## 2 1961 -0.234721502
## 3 1962 0.544102425
## 4 1963 -0.837212701
## 5 1964 -0.956356193
## 6 1965 -0.699724801
## 7 1966 -0.850601763
## 8 1967 -0.387628509
## 9 1968 0.061043154
## 10 1969 0.381984271
## 11 1970 -0.659303726
## 12 1971 0.126585564
## 13 1972 -0.009310217
## 14 1973 0.798500537
## 15 1974 1.863211350
## 16 1975 1.249270341
## 17 1976 0.504466636
## 18 1977 -0.737188575
## 19 1978 -0.749107276
## 20 1979 -0.221840562
## 21 1980 -0.686691381
## 22 1981 -0.092064694
## 23 1982 0.620750037
## 24 1983 -0.393901047
## 25 1984 -0.732847288
## 26 1985 0.429174729
## 27 1986 1.323636763
## 28 1987 0.612240617
## 29 1988 -0.201546121
## 30 1989 -0.498914711
## 31 1990 -1.094918870
## 32 1991 -1.901992588
## 33 1992 -1.211598783
## 34 1993 -1.039208546
## 35 1994 -0.972836612
## 36 1995 -0.657690616
## 37 1996 0.544553283
## 38 1997 1.515780517
## 39 1998 1.885861292
## 40 1999 2.080693375
## 41 2000 1.473593960
## 42 2001 0.958507189
## 43 2002 0.209460249
## 44 2003 -1.381917947
## 45 2004 -0.555041996
## 46 2005 0.593098207
## 47 2006 1.455274563
## 48 2007 0.736914838
## GR_IR GRCAT GRLOS GRMIN GRUMA GRWEN MFBEA MFBIG MFCAM
## [1,] 0.2112462 -0.03420949 0.1016297 0.1658748 0.1063399 0.07647642 0.3068064 0.3254104 0.3094061
## MFLOO MFMAR MFSUL SFEFS SFMAI SFSEC SREFS SRLEM SRLMA
## [1,] 0.3418419 0.1306568 0.2288056 0.2960428 0.2557321 0.4550384 0.3679295 0.1945687 0.3956831
## SRUMA SRVAL SRYFS
## [1,] 0.4925672 0.3156856 0.38948582.3 Risk-free index
For these purposes, the risk free index is zero for all years. The idea is that any population with a standardized return of zero is replacing itself and therefore not declining.
free <- rep(0, TT)3 CAPM in MARSS
3.1 v1: a static; b unconstrained
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## lower R block for beta
bb <- aa + nn
## equal var-cov
# QQ[bb,bb] <- "beta.cov"
# diag(QQ[bb,bb]) <- "beta.var"
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j+nn,k+nn] <- paste0("beta",j,k)
QQ[k+nn,j+nn] <- paste0("beta",j,k)
}
}
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1500 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1500 iterations.
## Log-likelihood: -390.0471
## AIC: 848.0942 AICc: 857.8085
##
## Estimate
## R.diag 0.706
## Q.beta1 0.719
## Q.beta12 0.942
## Q.beta13 0.756
## Q.beta14 0.604
## Q.beta15 0.872
## Q.beta16 0.762
## Q.beta2 1.234
## Q.beta23 0.991
## Q.beta24 0.791
## Q.beta25 1.142
## Q.beta26 0.999
## Q.beta3 0.796
## Q.beta34 0.635
## Q.beta35 0.917
## Q.beta36 0.802
## Q.beta4 0.507
## Q.beta45 0.732
## Q.beta46 0.640
## Q.beta5 1.057
## Q.beta56 0.924
## Q.beta6 0.808
## x0.X1 -0.287
## x0.X2 -0.326
## x0.X3 -0.158
## x0.X4 -0.174
## x0.X5 -0.344
## x0.X6 -0.293
## x0.X7 -0.743
## x0.X8 -1.343
## x0.X9 -0.880
## x0.X10 -0.614
## x0.X11 -1.034
## x0.X12 -0.943
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1374 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1374 iterations.
## Log-likelihood: -397.2663
## AIC: 862.5327 AICc: 872.7036
##
## Estimate
## R.diag 0.623
## Q.beta1 5.100
## Q.beta12 5.551
## Q.beta13 3.435
## Q.beta14 3.424
## Q.beta15 5.526
## Q.beta16 4.243
## Q.beta2 6.042
## Q.beta23 3.743
## Q.beta24 3.731
## Q.beta25 6.014
## Q.beta26 4.616
## Q.beta3 2.394
## Q.beta34 2.379
## Q.beta35 3.707
## Q.beta36 2.830
## Q.beta4 2.364
## Q.beta45 3.697
## Q.beta46 2.824
## Q.beta5 5.990
## Q.beta56 4.602
## Q.beta6 3.538
## x0.X1 -0.333
## x0.X2 -0.374
## x0.X3 -0.451
## x0.X4 -0.403
## x0.X5 -0.358
## x0.X6 -0.365
## x0.X7 -0.358
## x0.X8 -0.333
## x0.X9 -0.288
## x0.X10 -0.230
## x0.X11 -0.446
## x0.X12 -0.294
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1013 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1013 iterations.
## Log-likelihood: -167.7804
## AIC: 361.5608 AICc: 364.5444
##
## Estimate
## R.diag 0.370
## Q.beta1 1.405
## Q.beta12 1.076
## Q.beta13 1.112
## Q.beta2 0.825
## Q.beta23 0.852
## Q.beta3 0.881
## x0.X1 -0.277
## x0.X2 -0.237
## x0.X3 -0.143
## x0.X4 0.825
## x0.X5 0.636
## x0.X6 0.825
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1276 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1276 iterations.
## Log-likelihood: -371.8351
## AIC: 811.6702 AICc: 821.4645
##
## Estimate
## R.diag 0.381
## Q.beta1 9.543
## Q.beta12 7.379
## Q.beta13 6.356
## Q.beta14 7.864
## Q.beta15 8.931
## Q.beta16 12.867
## Q.beta2 5.796
## Q.beta23 4.951
## Q.beta24 6.078
## Q.beta25 6.838
## Q.beta26 9.838
## Q.beta3 4.249
## Q.beta34 5.238
## Q.beta35 5.923
## Q.beta36 8.528
## Q.beta4 6.482
## Q.beta45 7.364
## Q.beta46 10.611
## Q.beta5 8.413
## Q.beta56 12.131
## Q.beta6 17.495
## x0.X1 -0.645
## x0.X2 -0.696
## x0.X3 -0.566
## x0.X4 -0.535
## x0.X5 -0.576
## x0.X6 -0.779
## x0.X7 -5.837
## x0.X8 -4.579
## x0.X9 -3.728
## x0.X10 -4.512
## x0.X11 -5.483
## x0.X12 -7.824
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.1.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.2 v1: a static; b equalvarcov
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## lower R block for beta
bb <- aa + nn
## equal var-cov
QQ[bb,bb] <- "beta.cov"
diag(QQ[bb,bb]) <- "beta.var"
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1415 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1415 iterations.
## Log-likelihood: -392.2293
## AIC: 814.4587 AICc: 816.2769
##
## Estimate
## R.diag 0.720
## Q.beta.var 0.833
## Q.beta.cov 0.833
## x0.X1 -0.296
## x0.X2 -0.269
## x0.X3 -0.157
## x0.X4 -0.211
## x0.X5 -0.315
## x0.X6 -0.293
## x0.X7 -0.879
## x0.X8 -1.175
## x0.X9 -0.969
## x0.X10 -0.903
## x0.X11 -0.967
## x0.X12 -1.022
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1338 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1338 iterations.
## Log-likelihood: -399.7483
## AIC: 829.4966 AICc: 831.3938
##
## Estimate
## R.diag 0.7063
## Q.beta.var 3.6593
## Q.beta.cov 3.6593
## x0.X1 -0.3051
## x0.X2 -0.3234
## x0.X3 -0.4413
## x0.X4 -0.4151
## x0.X5 -0.3175
## x0.X6 -0.3518
## x0.X7 -0.1535
## x0.X8 -0.0502
## x0.X9 -0.1131
## x0.X10 -0.0547
## x0.X11 -0.1856
## x0.X12 -0.2453
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 984 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 984 iterations.
## Log-likelihood: -169.5842
## AIC: 357.1684 AICc: 358.597
##
## Estimate
## R.diag 0.397
## Q.beta.var 0.908
## Q.beta.cov 0.908
## x0.X1 -0.246
## x0.X2 -0.239
## x0.X3 -0.141
## x0.X4 0.712
## x0.X5 0.645
## x0.X6 0.818
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1458 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1458 iterations.
## Log-likelihood: -387.8977
## AIC: 805.7954 AICc: 807.6275
##
## Estimate
## R.diag 0.522
## Q.beta.var 6.799
## Q.beta.cov 6.799
## x0.X1 -0.505
## x0.X2 -0.611
## x0.X3 -0.600
## x0.X4 -0.507
## x0.X5 -0.549
## x0.X6 -0.635
## x0.X7 -3.115
## x0.X8 -3.437
## x0.X9 -3.266
## x0.X10 -3.025
## x0.X11 -3.189
## x0.X12 -2.944
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.2.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.3 v1: a static; b diag & unequal
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## lower R block for beta
bb <- aa + nn
## diag & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1534 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1534 iterations.
## Log-likelihood: -425.1244
## AIC: 888.2489 AICc: 891.172
##
## Estimate
## R.diag 1.11e+00
## Q.beta1 1.56e-02
## Q.beta2 8.14e-02
## Q.beta3 1.78e-02
## Q.beta4 4.08e-05
## Q.beta5 2.87e-02
## Q.beta6 1.08e-04
## x0.X1 -1.94e-01
## x0.X2 -2.68e-01
## x0.X3 -3.38e-02
## x0.X4 -2.97e-02
## x0.X5 -2.33e-01
## x0.X6 -1.18e-01
## x0.X7 -2.51e-01
## x0.X8 -1.19e+00
## x0.X9 -3.06e-02
## x0.X10 1.80e-01
## x0.X11 -7.18e-01
## x0.X12 7.81e-02
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 3014 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 3014 iterations.
## Log-likelihood: -452.207
## AIC: 942.414 AICc: 945.4662
##
## Estimate
## R.diag 1.65e+00
## Q.beta1 4.15e-05
## Q.beta2 3.08e-03
## Q.beta3 2.30e-02
## Q.beta4 1.46e-04
## Q.beta5 5.88e-05
## Q.beta6 4.99e-05
## x0.X1 -1.13e-01
## x0.X2 -1.49e-01
## x0.X3 -3.36e-01
## x0.X4 -1.21e-01
## x0.X5 -1.53e-01
## x0.X6 -3.63e-01
## x0.X7 4.13e-01
## x0.X8 3.58e-01
## x0.X9 -7.52e-02
## x0.X10 5.15e-01
## x0.X11 1.84e-01
## x0.X12 3.57e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 2082 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2082 iterations.
## Log-likelihood: -181.13
## AIC: 382.2601 AICc: 384.0201
##
## Estimate
## R.diag 0.817099
## Q.beta1 0.008137
## Q.beta2 0.000028
## Q.beta3 0.000066
## x0.X1 -0.166970
## x0.X2 -0.120705
## x0.X3 0.005672
## x0.X4 -0.048574
## x0.X5 0.229948
## x0.X6 0.485909
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1740 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1740 iterations.
## Log-likelihood: -463.9245
## AIC: 965.849 AICc: 968.7947
##
## Estimate
## R.diag 1.65e+00
## Q.beta1 1.95e-04
## Q.beta2 4.43e-05
## Q.beta3 4.26e-05
## Q.beta4 4.39e-05
## Q.beta5 7.75e-05
## Q.beta6 1.89e-04
## x0.X1 -1.43e-01
## x0.X2 -2.48e-01
## x0.X3 -1.86e-01
## x0.X4 -1.44e-01
## x0.X5 -2.04e-01
## x0.X6 -3.86e-01
## x0.X7 5.58e-01
## x0.X8 2.43e-01
## x0.X9 4.70e-01
## x0.X10 6.56e-01
## x0.X11 4.73e-01
## x0.X12 7.64e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.3.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.4 v2: a unconstrained; b unconstrained
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
# QQ[aa,aa] <- "alpha.cov"
# diag(QQ[aa,aa]) <- "alpha.var"
## diagonal & unequal
diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j,k] <- paste0("alpha",j,k)
QQ[k,j] <- paste0("alpha",j,k)
}
}
## lower R block for beta
bb <- aa + nn
## equal var-cov
# QQ[bb,bb] <- "beta.cov"
# diag(QQ[bb,bb]) <- "beta.var"
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j+nn,k+nn] <- paste0("beta",j,k)
QQ[k+nn,j+nn] <- paste0("beta",j,k)
}
}
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2703 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2703 iterations.
## Log-likelihood: -363.2693
## AIC: 836.5386 AICc: 864.0386
##
## Estimate
## R.diag 0.3003
## Q.alpha1 0.4947
## Q.alpha12 0.5072
## Q.alpha13 0.6028
## Q.alpha14 0.4998
## Q.alpha15 0.2693
## Q.alpha16 0.4113
## Q.alpha2 0.7917
## Q.alpha23 0.5354
## Q.alpha24 0.3092
## Q.alpha25 0.3237
## Q.alpha26 0.3896
## Q.alpha3 0.7968
## Q.alpha34 0.6822
## Q.alpha35 0.4049
## Q.alpha36 0.5219
## Q.alpha4 0.6614
## Q.alpha45 0.2649
## Q.alpha46 0.4461
## Q.alpha5 0.3794
## Q.alpha56 0.2462
## Q.alpha6 0.3588
## Q.beta1 0.0132
## Q.beta12 0.0388
## Q.beta13 0.0140
## Q.beta14 0.0324
## Q.beta15 0.0979
## Q.beta16 0.0535
## Q.beta2 0.1244
## Q.beta23 0.0465
## Q.beta24 0.1163
## Q.beta25 0.3303
## Q.beta26 0.1848
## Q.beta3 0.0177
## Q.beta34 0.0454
## Q.beta35 0.1260
## Q.beta36 0.0711
## Q.beta4 0.1226
## Q.beta45 0.3271
## Q.beta46 0.1876
## Q.beta5 0.9014
## Q.beta56 0.5103
## Q.beta6 0.2903
## x0.X1 -0.4679
## x0.X2 -1.7161
## x0.X3 0.3303
## x0.X4 0.5274
## x0.X5 0.0813
## x0.X6 -0.5697
## x0.X7 -0.1207
## x0.X8 -0.9706
## x0.X9 -0.0274
## x0.X10 -0.0442
## x0.X11 -1.5186
## x0.X12 -0.6445
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2135 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2135 iterations.
## Log-likelihood: -352.7901
## AIC: 815.5802 AICc: 844.5004
##
## Estimate
## R.diag 0.281363
## Q.alpha1 1.429407
## Q.alpha12 1.613941
## Q.alpha13 0.673732
## Q.alpha14 0.692343
## Q.alpha15 1.489888
## Q.alpha16 1.608054
## Q.alpha2 1.842574
## Q.alpha23 0.797055
## Q.alpha24 0.823314
## Q.alpha25 1.681005
## Q.alpha26 1.776505
## Q.alpha3 0.924359
## Q.alpha34 1.097678
## Q.alpha35 0.720479
## Q.alpha36 0.109656
## Q.alpha4 1.317354
## Q.alpha45 0.745323
## Q.alpha46 -0.045377
## Q.alpha5 1.553820
## Q.alpha56 1.656737
## Q.alpha6 2.501869
## Q.beta1 0.002804
## Q.beta12 0.001372
## Q.beta13 -0.002389
## Q.beta14 0.002082
## Q.beta15 0.006010
## Q.beta16 -0.006659
## Q.beta2 0.000708
## Q.beta23 -0.001106
## Q.beta24 0.001024
## Q.beta25 0.002905
## Q.beta26 -0.003139
## Q.beta3 0.002238
## Q.beta34 -0.001723
## Q.beta35 -0.005229
## Q.beta36 0.005939
## Q.beta4 0.001635
## Q.beta45 0.004485
## Q.beta46 -0.005042
## Q.beta5 0.012994
## Q.beta56 -0.014538
## Q.beta6 0.016617
## x0.X1 0.231601
## x0.X2 -0.288871
## x0.X3 0.008899
## x0.X4 0.223452
## x0.X5 0.335920
## x0.X6 0.100134
## x0.X7 0.509994
## x0.X8 0.522791
## x0.X9 0.094360
## x0.X10 0.126260
## x0.X11 0.478596
## x0.X12 0.416539
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 2136 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2136 iterations.
## Log-likelihood: -135.2344
## AIC: 308.4688 AICc: 315.0206
##
## Estimate
## R.diag 1.37e-01
## Q.alpha1 7.71e-01
## Q.alpha12 6.34e-01
## Q.alpha13 6.02e-01
## Q.alpha2 5.23e-01
## Q.alpha23 4.89e-01
## Q.alpha3 5.48e-01
## Q.beta1 1.80e-04
## Q.beta12 1.01e-04
## Q.beta13 1.06e-04
## Q.beta2 6.48e-05
## Q.beta23 6.18e-05
## Q.beta3 7.57e-05
## x0.X1 -9.75e-01
## x0.X2 -8.61e-01
## x0.X3 -6.27e-01
## x0.X4 3.00e-01
## x0.X5 2.41e-01
## x0.X6 4.25e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1577 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1577 iterations.
## Log-likelihood: -338.4211
## AIC: 786.8422 AICc: 814.59
##
## Estimate
## R.diag 0.28630
## Q.alpha1 1.34857
## Q.alpha12 0.96660
## Q.alpha13 0.95088
## Q.alpha14 1.06089
## Q.alpha15 1.28876
## Q.alpha16 1.68622
## Q.alpha2 0.69572
## Q.alpha23 0.67684
## Q.alpha24 0.77004
## Q.alpha25 0.93189
## Q.alpha26 1.22611
## Q.alpha3 0.67868
## Q.alpha34 0.73401
## Q.alpha35 0.89737
## Q.alpha36 1.16775
## Q.alpha4 0.87348
## Q.alpha45 1.04886
## Q.alpha46 1.41327
## Q.alpha5 1.26369
## Q.alpha56 1.69375
## Q.alpha6 2.33465
## Q.beta1 0.00653
## Q.beta12 0.03641
## Q.beta13 0.01717
## Q.beta14 0.00838
## Q.beta15 -0.01926
## Q.beta16 -0.01446
## Q.beta2 0.20737
## Q.beta23 0.09574
## Q.beta24 0.04554
## Q.beta25 -0.11392
## Q.beta26 -0.08513
## Q.beta3 0.04585
## Q.beta34 0.02275
## Q.beta35 -0.04938
## Q.beta36 -0.03737
## Q.beta4 0.01184
## Q.beta45 -0.02161
## Q.beta46 -0.01663
## Q.beta5 0.06896
## Q.beta56 0.05064
## Q.beta6 0.03738
## x0.X1 -0.08959
## x0.X2 -0.17092
## x0.X3 -0.27497
## x0.X4 0.36118
## x0.X5 0.38791
## x0.X6 1.08623
## x0.X7 0.49013
## x0.X8 0.10929
## x0.X9 0.49369
## x0.X10 0.71621
## x0.X11 0.66461
## x0.X12 0.88863
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.4.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.5 v2: a unconstrained; b equalvarcov
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
# QQ[aa,aa] <- "alpha.cov"
# diag(QQ[aa,aa]) <- "alpha.var"
## diagonal & unequal
diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j,k] <- paste0("alpha",j,k)
QQ[k,j] <- paste0("alpha",j,k)
}
}
## lower R block for beta
bb <- aa + nn
## equal var-cov
QQ[bb,bb] <- "beta.cov"
diag(QQ[bb,bb]) <- "beta.var"
## diagonal & unequal
# diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
# for(j in 1:(nn-1)) {
# for(k in (j+1):nn) {
# QQ[j+nn,k+nn] <- paste0("beta",j,k)
# QQ[k+nn,j+nn] <- paste0("beta",j,k)
# }
# }
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1956 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1956 iterations.
## Log-likelihood: -367.3308
## AIC: 806.6615 AICc: 817.6245
##
## Estimate
## R.diag 4.70e-01
## Q.alpha1 4.69e-01
## Q.alpha12 5.73e-01
## Q.alpha13 5.86e-01
## Q.alpha14 4.67e-01
## Q.alpha15 5.27e-01
## Q.alpha16 4.84e-01
## Q.alpha2 8.52e-01
## Q.alpha23 6.45e-01
## Q.alpha24 4.34e-01
## Q.alpha25 6.06e-01
## Q.alpha26 5.36e-01
## Q.alpha3 7.73e-01
## Q.alpha34 6.50e-01
## Q.alpha35 6.87e-01
## Q.alpha36 6.25e-01
## Q.alpha4 5.89e-01
## Q.alpha45 5.63e-01
## Q.alpha46 5.29e-01
## Q.alpha5 6.15e-01
## Q.alpha56 5.51e-01
## Q.alpha6 5.20e-01
## Q.beta.var 6.89e-05
## Q.beta.cov 5.03e-05
## x0.X1 -3.95e-01
## x0.X2 -1.24e+00
## x0.X3 2.67e-01
## x0.X4 5.05e-01
## x0.X5 2.21e-02
## x0.X6 -2.13e-01
## x0.X7 2.70e-01
## x0.X8 -1.02e-01
## x0.X9 2.08e-01
## x0.X10 2.88e-01
## x0.X11 1.81e-01
## x0.X12 1.64e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2051 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2051 iterations.
## Log-likelihood: -353.0967
## AIC: 778.1934 AICc: 789.6761
##
## Estimate
## R.diag 3.06e-01
## Q.alpha1 1.44e+00
## Q.alpha12 1.59e+00
## Q.alpha13 6.86e-01
## Q.alpha14 7.06e-01
## Q.alpha15 1.50e+00
## Q.alpha16 1.63e+00
## Q.alpha2 1.78e+00
## Q.alpha23 8.02e-01
## Q.alpha24 8.34e-01
## Q.alpha25 1.66e+00
## Q.alpha26 1.76e+00
## Q.alpha3 9.56e-01
## Q.alpha34 1.11e+00
## Q.alpha35 7.30e-01
## Q.alpha36 1.43e-01
## Q.alpha4 1.29e+00
## Q.alpha45 7.54e-01
## Q.alpha46 2.23e-02
## Q.alpha5 1.57e+00
## Q.alpha56 1.69e+00
## Q.alpha6 2.49e+00
## Q.beta.var 6.00e-05
## Q.beta.cov 4.71e-05
## x0.X1 1.94e-01
## x0.X2 -2.56e-01
## x0.X3 1.09e-01
## x0.X4 1.72e-01
## x0.X5 2.56e-01
## x0.X6 1.56e-01
## x0.X7 4.48e-01
## x0.X8 4.99e-01
## x0.X9 1.64e-01
## x0.X10 1.01e-01
## x0.X11 3.59e-01
## x0.X12 5.38e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1564 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1564 iterations.
## Log-likelihood: -135.2571
## AIC: 300.5141 AICc: 304.5141
##
## Estimate
## R.diag 1.37e-01
## Q.alpha1 7.77e-01
## Q.alpha12 6.35e-01
## Q.alpha13 6.04e-01
## Q.alpha2 5.19e-01
## Q.alpha23 4.87e-01
## Q.alpha3 5.48e-01
## Q.beta.var 4.78e-05
## Q.beta.cov 3.70e-05
## x0.X1 -9.73e-01
## x0.X2 -8.56e-01
## x0.X3 -6.27e-01
## x0.X4 3.09e-01
## x0.X5 2.44e-01
## x0.X6 4.31e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 2028 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2028 iterations.
## Log-likelihood: -340.9977
## AIC: 753.9953 AICc: 765.0492
##
## Estimate
## R.diag 3.56e-01
## Q.alpha1 1.32e+00
## Q.alpha12 1.00e+00
## Q.alpha13 9.35e-01
## Q.alpha14 1.07e+00
## Q.alpha15 1.23e+00
## Q.alpha16 1.67e+00
## Q.alpha2 7.79e-01
## Q.alpha23 7.13e-01
## Q.alpha24 7.91e-01
## Q.alpha25 8.92e-01
## Q.alpha26 1.20e+00
## Q.alpha3 6.63e-01
## Q.alpha34 7.52e-01
## Q.alpha35 8.63e-01
## Q.alpha36 1.17e+00
## Q.alpha4 8.80e-01
## Q.alpha45 1.03e+00
## Q.alpha46 1.42e+00
## Q.alpha5 1.23e+00
## Q.alpha56 1.71e+00
## Q.alpha6 2.39e+00
## Q.beta.var 8.35e-05
## Q.beta.cov 6.97e-05
## x0.X1 6.78e-03
## x0.X2 -4.22e-01
## x0.X3 -1.70e-01
## x0.X4 2.46e-01
## x0.X5 5.47e-01
## x0.X6 9.87e-01
## x0.X7 5.06e-01
## x0.X8 1.83e-01
## x0.X9 3.89e-01
## x0.X10 6.20e-01
## x0.X11 4.60e-01
## x0.X12 6.94e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.5.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.6 v3: a equalvarcov; b unconstrained
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
QQ[aa,aa] <- "alpha.cov"
diag(QQ[aa,aa]) <- "alpha.var"
## lower R block for beta
bb <- aa + nn
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j+nn,k+nn] <- paste0("beta",j,k)
QQ[k+nn,j+nn] <- paste0("beta",j,k)
}
}
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1962 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1962 iterations.
## Log-likelihood: -371.5489
## AIC: 815.0978 AICc: 826.0608
##
## Estimate
## R.diag 0.512873
## Q.alpha.var 0.472812
## Q.alpha.cov 0.472788
## Q.beta1 0.012384
## Q.beta12 0.033173
## Q.beta13 -0.000533
## Q.beta14 -0.014850
## Q.beta15 0.005518
## Q.beta16 -0.001409
## Q.beta2 0.089998
## Q.beta23 -0.003417
## Q.beta24 -0.048698
## Q.beta25 0.001972
## Q.beta26 -0.012209
## Q.beta3 0.004044
## Q.beta34 0.018773
## Q.beta35 0.026132
## Q.beta36 0.017356
## Q.beta4 0.100114
## Q.beta45 0.113228
## Q.beta46 0.080249
## Q.beta5 0.177245
## Q.beta56 0.113874
## Q.beta6 0.075208
## x0.X1 -0.271026
## x0.X2 -0.350586
## x0.X3 -0.071601
## x0.X4 -0.085512
## x0.X5 -0.338093
## x0.X6 -0.222222
## x0.X7 -0.158583
## x0.X8 -1.184148
## x0.X9 0.055420
## x0.X10 0.015174
## x0.X11 -1.186145
## x0.X12 -0.568935
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 4799 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 4799 iterations.
## Log-likelihood: -376.5283
## AIC: 825.0567 AICc: 836.5394
##
## Estimate
## R.diag 0.574195
## Q.alpha.var 1.090505
## Q.alpha.cov 1.090496
## Q.beta1 0.001277
## Q.beta12 0.002277
## Q.beta13 0.010294
## Q.beta14 0.009033
## Q.beta15 -0.000374
## Q.beta16 -0.005056
## Q.beta2 0.004094
## Q.beta23 0.018539
## Q.beta24 0.016265
## Q.beta25 -0.000690
## Q.beta26 -0.009128
## Q.beta3 0.084323
## Q.beta34 0.073969
## Q.beta35 -0.003216
## Q.beta36 -0.041622
## Q.beta4 0.064896
## Q.beta45 -0.002819
## Q.beta46 -0.036510
## Q.beta5 0.000152
## Q.beta56 0.001618
## Q.beta6 0.020596
## x0.X1 0.159421
## x0.X2 0.136142
## x0.X3 -0.000687
## x0.X4 0.044807
## x0.X5 0.158152
## x0.X6 0.105142
## x0.X7 0.433713
## x0.X8 0.531834
## x0.X9 0.377478
## x0.X10 0.415484
## x0.X11 0.342442
## x0.X12 0.325962
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 2022 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2022 iterations.
## Log-likelihood: -139.7712
## AIC: 309.5423 AICc: 313.5423
##
## Estimate
## R.diag 0.20940
## Q.alpha.var 0.50907
## Q.alpha.cov 0.50906
## Q.beta1 0.00870
## Q.beta12 0.00481
## Q.beta13 0.00492
## Q.beta2 0.00267
## Q.beta23 0.00272
## Q.beta3 0.00279
## x0.X1 -0.90546
## x0.X2 -0.91729
## x0.X3 -0.78593
## x0.X4 -0.06701
## x0.X5 0.06959
## x0.X6 0.22662
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1864 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1864 iterations.
## Log-likelihood: -352.8772
## AIC: 777.7545 AICc: 788.8084
##
## Estimate
## R.diag 0.38467
## Q.alpha.var 0.97030
## Q.alpha.cov 0.97029
## Q.beta1 0.00124
## Q.beta12 -0.00892
## Q.beta13 -0.00240
## Q.beta14 0.00136
## Q.beta15 0.01100
## Q.beta16 0.02011
## Q.beta2 0.08153
## Q.beta23 0.01662
## Q.beta24 -0.02014
## Q.beta25 -0.11016
## Q.beta26 -0.20117
## Q.beta3 0.00545
## Q.beta34 -0.00109
## Q.beta35 -0.01852
## Q.beta36 -0.03411
## Q.beta4 0.00940
## Q.beta45 0.03298
## Q.beta46 0.05982
## Q.beta5 0.15644
## Q.beta56 0.28511
## Q.beta6 0.51972
## x0.X1 0.18842
## x0.X2 0.05143
## x0.X3 0.12064
## x0.X4 0.22757
## x0.X5 0.24395
## x0.X6 0.17975
## x0.X7 0.42536
## x0.X8 0.46287
## x0.X9 0.54532
## x0.X10 0.69797
## x0.X11 0.24589
## x0.X12 0.36009
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.6.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.7 v3: a equalvarcov; b equalvarcov
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
QQ[aa,aa] <- "alpha.cov"
diag(QQ[aa,aa]) <- "alpha.var"
## lower R block for beta
bb <- aa + nn
## equal var-cov
QQ[bb,bb] <- "beta.cov"
diag(QQ[bb,bb]) <- "beta.var"
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2023 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2023 iterations.
## Log-likelihood: -371.9797
## AIC: 777.9595 AICc: 780.2954
##
## Estimate
## R.diag 6.00e-01
## Q.alpha.var 5.16e-01
## Q.alpha.cov 5.16e-01
## Q.beta.var 6.42e-05
## Q.beta.cov 4.46e-05
## x0.X1 -1.94e-01
## x0.X2 -1.64e-01
## x0.X3 -5.67e-02
## x0.X4 -1.12e-01
## x0.X5 -2.22e-01
## x0.X6 -1.76e-01
## x0.X7 2.85e-01
## x0.X8 -7.16e-03
## x0.X9 1.94e-01
## x0.X10 2.58e-01
## x0.X11 1.92e-01
## x0.X12 1.46e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1943 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1943 iterations.
## Log-likelihood: -378.4222
## AIC: 790.8445 AICc: 793.2827
##
## Estimate
## R.diag 6.29e-01
## Q.alpha.var 1.12e+00
## Q.alpha.cov 1.12e+00
## Q.beta.var 7.72e-05
## Q.beta.cov 5.71e-05
## x0.X1 1.42e-01
## x0.X2 1.23e-01
## x0.X3 1.12e-02
## x0.X4 6.52e-02
## x0.X5 1.45e-01
## x0.X6 9.14e-02
## x0.X7 3.97e-01
## x0.X8 5.01e-01
## x0.X9 4.26e-01
## x0.X10 4.81e-01
## x0.X11 3.22e-01
## x0.X12 2.90e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1688 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1688 iterations.
## Log-likelihood: -140.0225
## AIC: 302.0449 AICc: 304.1739
##
## Estimate
## R.diag 2.15e-01
## Q.alpha.var 5.30e-01
## Q.alpha.cov 5.30e-01
## Q.beta.var 7.23e-05
## Q.beta.cov 5.88e-05
## x0.X1 -8.70e-01
## x0.X2 -8.96e-01
## x0.X3 -7.63e-01
## x0.X4 3.27e-01
## x0.X5 2.65e-01
## x0.X6 4.36e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1995 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1995 iterations.
## Log-likelihood: -358.7276
## AIC: 751.4551 AICc: 753.809
##
## Estimate
## R.diag 4.92e-01
## Q.alpha.var 1.00e+00
## Q.alpha.cov 1.00e+00
## Q.beta.var 6.84e-05
## Q.beta.cov 5.03e-05
## x0.X1 2.22e-01
## x0.X2 1.16e-01
## x0.X3 1.29e-01
## x0.X4 2.23e-01
## x0.X5 1.86e-01
## x0.X6 9.38e-02
## x0.X7 5.22e-01
## x0.X8 1.97e-01
## x0.X9 3.74e-01
## x0.X10 6.10e-01
## x0.X11 4.37e-01
## x0.X12 6.92e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.7.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.8 v4: a diag & unequal; b diag & unequal
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
# QQ[aa,aa] <- "alpha.cov"
# diag(QQ[aa,aa]) <- "alpha.var"
## diagonal & unequal
diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
## lower R block for beta
bb <- aa + nn
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2398 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2398 iterations.
## Log-likelihood: -423.8285
## AIC: 897.657 AICc: 902.7751
##
## Estimate
## R.diag 9.57e-01
## Q.alpha1 2.50e-05
## Q.alpha2 3.19e-01
## Q.alpha3 1.63e-01
## Q.alpha4 4.52e-05
## Q.alpha5 2.91e-05
## Q.alpha6 1.94e-02
## Q.beta1 1.81e-02
## Q.beta2 1.95e-04
## Q.beta3 4.10e-05
## Q.beta4 2.18e-05
## Q.beta5 3.96e-02
## Q.beta6 3.00e-05
## x0.X1 -2.03e-01
## x0.X2 -8.69e-01
## x0.X3 7.74e-01
## x0.X4 -3.23e-02
## x0.X5 -2.51e-01
## x0.X6 -2.97e-01
## x0.X7 -2.82e-01
## x0.X8 -1.24e-01
## x0.X9 1.22e-01
## x0.X10 1.80e-01
## x0.X11 -8.41e-01
## x0.X12 6.52e-02
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1998 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1998 iterations.
## Log-likelihood: -452.9958
## AIC: 955.9916 AICc: 961.3414
##
## Estimate
## R.diag 1.54e+00
## Q.alpha1 7.87e-05
## Q.alpha2 2.98e-01
## Q.alpha3 5.50e-05
## Q.alpha4 1.44e-04
## Q.alpha5 6.11e-05
## Q.alpha6 6.17e-05
## Q.beta1 7.04e-05
## Q.beta2 8.90e-05
## Q.beta3 2.55e-02
## Q.beta4 3.54e-03
## Q.beta5 9.93e-05
## Q.beta6 8.15e-05
## x0.X1 -1.17e-01
## x0.X2 -8.48e-01
## x0.X3 -3.39e-01
## x0.X4 -1.57e-01
## x0.X5 -1.55e-01
## x0.X6 -3.63e-01
## x0.X7 4.09e-01
## x0.X8 4.56e-01
## x0.X9 -8.80e-02
## x0.X10 3.41e-01
## x0.X11 1.81e-01
## x0.X12 3.54e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1257 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1257 iterations.
## Log-likelihood: -177.6811
## AIC: 381.3621 AICc: 384.3457
##
## Estimate
## R.diag 2.24e-01
## Q.alpha1 6.78e-01
## Q.alpha2 2.49e-01
## Q.alpha3 3.72e-01
## Q.beta1 1.58e-04
## Q.beta2 5.87e-05
## Q.beta3 1.02e-04
## x0.X1 -3.95e-01
## x0.X2 -1.22e+00
## x0.X3 -5.32e-01
## x0.X4 2.65e-01
## x0.X5 2.16e-01
## x0.X6 4.79e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1405 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1405 iterations.
## Log-likelihood: -458.2313
## AIC: 966.4625 AICc: 971.6213
##
## Estimate
## R.diag 5.73e-01
## Q.alpha1 7.26e-01
## Q.alpha2 4.67e-01
## Q.alpha3 2.71e-01
## Q.alpha4 3.54e-01
## Q.alpha5 8.82e-01
## Q.alpha6 1.66e+00
## Q.beta1 8.31e-05
## Q.beta2 6.73e-05
## Q.beta3 6.61e-05
## Q.beta4 5.74e-05
## Q.beta5 1.02e-04
## Q.beta6 2.26e-04
## x0.X1 -1.05e-01
## x0.X2 -4.65e-01
## x0.X3 -4.28e-01
## x0.X4 2.56e-01
## x0.X5 1.23e-01
## x0.X6 1.03e+00
## x0.X7 5.02e-01
## x0.X8 1.90e-01
## x0.X9 4.74e-01
## x0.X10 6.16e-01
## x0.X11 4.48e-01
## x0.X12 9.61e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.8.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)